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ABSTRACT 

Motivated by path-integral numerical solutions of diffusion processes, PATHINT, we present 
a new tree algorithm, PATHTREE, which permits extremely fast accurate computation of 
probability distributions of a large class of general nonlinear diffusion processes. 
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1. INTRODUCTION 



1.1. Path Integral Solution of Diffusion Processes 

There are three equivalent mathematical representations of a diffusion process, provided of course 
that boundary and initial conditions can be properly specified in each representation. In this paper we 
refer to all three representations. 

The Langevin rate equation for a stochastic process in dS can be written as in a prepoint 
discretization, 

dS = fdt + gdW , 



<dW>=0, 



< dW{t)dW{f) >= dtS{t -f), (1) 

for general drift / and standard deviation g which may depend on S and t, wherein / and g are 
understood to be evaluated at the prepoint t. Here, we just consider S dependent, but our algorithm can 
easily be extended to time dependent cases and to multivariate systems. 

This corresponds to a Fokker-Planck equation representing the short-time conditional probability P 
of evolving within time dt, 

dP ^ d(JP) 1 d\g^P) 

dt 35 2 352 ' ^ ^ 

where the diffusion is given by g . 

The path-integral representation for P for the short-time propagator is given by 

P{S', t'lS, t) = &xp{-Ut) 



2g2 
dS S'-S 



,dt = t'-t. (3) 



dt dt 

In the above we have expUcitly used the prepoint discretization [1]. 



1.2. PATHINT Motivation From Previous Study 

In the above we have explicitly used the prepoint discretization, wherein / and g are understood to 
be evaluated at the prepoint t. In this paper, we do not require multivariate generalizations, or issues 
dealing with other discretizations, or explication of long-time path-integral evaluations, or issues dealing 
with Riemannian invariance of our distributions. There exist other references deaUng with these issues in 
the context of calculations presented here [2-5]. 

Our approach is motivated by a multivariable generahzation of a numerical path-integral 
algorithm [6-8], PATHINT, used to develop the long-time evolution of the short-time probability 
distribution as used in several studies in chaotic systems [9,10], neuroscience [9,11,12], and financial 
markets [4]. These studies suggested that we apply some aspects of this algorithm to the standard 
binomial tree. 



1.3. PATHTREE Algorithms 

Tree algorithms are generally derived from binomial random walks [13]. For many applications, 
"tree" algorithms are often used, corresponding to the above Langevin and Fokker-Planck 
equations [14,15]. These algorithms have typically been only well defined for specific functional forms of 
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faadg. 

We have previously presented a powerful PATHINT algorithm to deal with quite general / and g 
functions [4]. This general PATHTREE algorithm can be used beyond previous specific systems, 
affording fast reasonable resolution calculations for probability distributions of a large class of nonlinear 
diffusion problems. 

1.4. Organization of Paper 

Section 2 describes the standard tree algorithm. Section 3 develops our probabiUty PATHTREE 
algorithm. Section 4 presents our probabiUty calculations. Section 5 is our conclusion. 

2. STANDARD OPTION TREE ALGORITHM 

2.1. Binomial Tree 

In a two-step binomial tree, the step up Su or step down Sd from a given node at S is chosen to 
match the standard deviation of the differential process. The constraints on u and d are chosen as 

ud = l, (4) 

If we assign probability p to the up step Su, and q = (I- p)to the down step Sd, the matched mean and 
variance are 

pSu + (1 - p)Sd =< S{t + At) > , 

S^ipu^ + qd^ - (pu + qdf) = < iS(t + At)-< S(t + At) >f > . (5) 
The right-hand-side can be derived from the stochastic model used. 

2.2. Trinomial Tree 

The trinomial tree can be used as a robust alternate to the binomial tree. Assume pj,, p„, and p^ are 
the probabilities of up jump Su, middle (no-)jump S and down jump Sd, where the jumps are chosen to 
match the standard deviation. To match the variance of the process, the equations must satisfy 

Pu + Pm + Pd = '^ , 

S{PuU + Pm + Pdd) =<S{t + At)>, 

S^ip^u^ + Pm + Pdd^ - (PuU + Pm + Pddf) = < iS(t + At)-< S(t + At) >)2 > . (6) 

3. PROBABILITY TREE ALGORITHM 

3.1. General Diffusion Process 

Consider the general Markov multiplicative diffusion process interpreted as an Ito prepoint 
discretized process, Eq. (1) with drift / and diffusion g^. For financial option studies the particular form 
of the drift bS and diffusion (aS) is chosen for lognormal Black-Scholes (BS) calculations [14]. For 
options, the coefficient b is the cost of carry, e.g., b = r, the risk-free rate, when 5 is a stock price, and 
b = when 5 is a futures price [16]. The case of drift bS and constant diffusion diffusion cr^ corresponds 
to the Omstein-Uhlenbeck (OU) process [17]. 

Our formalism is general and can be applied to other functional forms of interest with quite general 
nonlinear drifts and diffusions, of course provided that the region of solution does not violate any 
boundary or initial conditions. 

Statistical properties of the dS process and of any derivative one based on nonlinear transformations 
apphed to S are determined once the transition probabiUty distribution function P(S, tlSo, to) is known, 
where the index denotes initial values of time and of the stochastic variable 5. Transformation are 
common and convenient for BS, 
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z = ln5, (7) 

yielding a simple Gaussian distribution in z, greatly simplifying analytic and numerical calculations. 

The probability distribution can be obtained by solving the associated forward Fokker-Planck 
equation Eq. (2). Appropriate boundaries and initial condition must be specified, e.g., 
PiS\So) = SiS-So). 

In general cases, the Fokker-Planck equation is rather difficult to solve, although a vast body of 
work is devoted to it [17]. The particular BS and OU cases possess exact results. 

Our goal is to obtain the solution of Eq. (1) for the more general process. A quite general code, 
PATHINT [4], works fine, but it is much slower than the PATHTREE algorithm we present here. 

3.2. Deficiency of Standard Algorithm to Order 'idt 

We briefly describe the CRR construction of the binomial tree approximation [18]. 

A tree is constructed that represents the time evolution of the stochastic variable S. 5 is assumed to 
take only 2 values, u, (up value), and d (down value) at moment t, given the value S at moment t — At. 
The probabilities for the up and down movements are p and q, respectively. The 4 unknowns {u, d, p, q} 
are calculated by imposing the normalization of the probability and matching the first two moments 
conditioned by the value S at t — At, using the variance of the exact probability distribution P(S, t\S(), to). 
One additional condition is arbitrary and is usually used to symmetrize the tree, e.g., ud = I. 

The main problem is that the above procedure cannot be applied to a general nonlinear diffusion 
process as considered in Eq. (1), as the algorithm involves a previous knowledge of terms of 0{At) in the 
formulas of quantities {u, p} obtained from a finite time expansion of the exact solution P sought. 
Otherwise, the discrete numerical approximation obtained does not converge to the proper solution. 

This observation can be checked analytically in the BS CRR case by replacing the relation 
u = exp(crAt) [15] with u = 1 + crfKt, and deriving the continuous limit of the tree. This also can be 
checked numerically, as when {u, p} are expanded to 0(At), the proper solution is obtained. 



3.3. Probability PATHTREE 

As mentioned previously, a general path-integral solution of the Fokker-Plank equation, including 
the Black-Scholes equation, can be numerically calculated using the PATHINT algorithm. Although this 
approach leads to very accurate results, it is computationally intensive. 

In order to obtain tree variables valid up to 0(At), we turn to the short-time path-integral 
representation of the solution of the Fokker-Planck equation, which is just the multiplicative Gaussian- 
Markovian distribution [1,19]. In the prepoint discretization relevant to the construction of a tree. 



P(S\ ns, t) = ^ exp 

^iKAtg^ 



(5' - 5 - fdtf 
Ig^At 



At = t'-t (8) 



valid for displacements 5" from S "reasonable" as measured by the standard deviation gVAf, which is the 
basis for the construction of meshes in the PATHINT algorithm. 

The crucial aspects of this approach are: There is no a priori need of the first moments of the exact 
long-time probabihty distribution P, as the necessary statistical information to the correct order in time is 
contained in the short-time propagator. The mesh in S at every time step need not recombine in the sense 
that the prepoint-postpoint relationship be the same among neighboring 5 nodes, as the short-time 
probabihty density gives the correct result up to order 0{At) for any final point S'. Instead, we use the 
natural metric of the space to first lay down our mesh, then dynamically calculate the evolving local short- 
time distributions on this mesh. 

We construct an additive PATHTREE, starting with the initial value 5o, with successive increments 
SM = Si + g<At ,Si>SQ 
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Si.i = Si-g<Ai ,Si<So, (9) 

where g is evaluated at 5,. We define the up and down probabilities p and q, resp., in an abbreviated 
notation, as 

P(i + l\i;At) 



^ P(i+l\i;At) + P(i-l\i;At) 
q=\- p . 



(10) 



where the short- time transition probabihty densities P's are calculated from Eq. (8). Note that in the limit 
of small At, 

lim;? = -. (11) 

Af->0 2 



3.3.1. Continuous Limit of the PATHTREE Algorithm 

In either the upper or lower branches of the tree (5, > or 5, < Sq, resp.), we always have the 
postpoint Sj+i in terms of the prepoint 5,, but we also need the inverses to find the asymptotic lim of our 

A->0 

PATHTREE building technique. For example, for the upper branch case. 



(12) 



This expression must be used to extract the down change d (5,_i = 5, + d), for comparison to the standard 
tree algorithm. 

The continuous limit of the previous tree building procedure is obtained by Taylor expanding all 
factors up to terms 0{At) and as functions of the prepoint 5, [15]. This leads to 

pu + qd 



At 



»f + 0(At"^) 



P'^'lf ^g^^OjAt^'^), 



(13) 



from which the correct partial differential equation, Eq. (2), is recovered up to 0{At). 

In implementing the PATHTREE algorithm, good numerical results are obtained in the parameter 
region defined by the convergence condition 



8 



9g 



dt + gidt 



'Si«\ 



(14) 



This insures the proper construction of the tree to order 0{At). 



3.3.2. Treating Binomial Oscillations 

Binomial trees exhibit by construction a systematic oscillatory behavior as a function of the number 
of steps in the tree (equivalently, the number of time slices used to build the tree), and the new building 
algorithm based on the short-time propagator of the path-integral representation of the solution of the 
Fokker-Planck equation has this same problem. A common practice [20] is to perform averages of runs 
with consecutive numbers of steps, e.g.. 



C = 



(15) 



where signifies the value calculated with number of steps. 
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3.3.3. Inappropriate TVinomial Tree 

Another type of tree is the trinomial tree discussed above, equivalent to the explicit finite difference 
method [14,15]. If we were to apply this approach to our new PATHTREE algorithm, we would allow the 
stochastic variable 5 to remain unchanged after time step A? with a certain probabiUty p^. However, in 
our construction the trinomial tree approach is not correct, as the deterministic path is dominant in the 
construction of the probabilities p^, p^}, and we would obtain 



3.4. Linear and Quadratic Aspects of Numerical Implementation 

PATHTREE computes the expected value of a random variable at a later time given the diffusion 
process and an initial point. The algorithm is similar to a binomial tree computation and it consists of 
three main steps: computation of a mesh of points, computation of transition probabiHties at those points 
and computation of the expected value of the random variable. 

The first step is the creation of a one dimensional mesh of points with gaps determined by the 
second moment of the short term distribution of the process. The mesh is created sequentially, starting 
from the initial point, by progressively adding to the last point already determined (for the upward part of 
the mesh) the value of the standard deviation of the short term distribution with the same point as 
prepoint. In a similar fashion we create the mesh downwards, this time by subtracting the standard 
deviations. The procedure takes a linear amount of time on the number of time slices being considered 
and contributes very little to the overall time of the algorithm. 

In the second step an array of up and down probabiUties is created. These probabilities are the 
values of the short term transition probability density function obtained by using the current point as 
prepoint and the two neighboring points as post points. The probabilities are renormalized to sum to 
unity. This procedure takes a linear amount of time on the number of time slices. Notice that the 
probabiUties only depend on the current point and not on time slice, hence only two probabilities are 
computed per element of the array of points. 

The third step is the computation of the expected value of the random variable. For example, the 
option price C is developed by marching backwards along the tree and applying the risk-neutral 
evaluation 



We emphasize again that in Ito terms the prepoint value is S;. This part works as a normal binomial tree 
algorithm. The algorithm uses the expected values at one time slice to compute the expected values at the 
previous one. The bulk of the time is spent in this part of the algorithm because the number of iterations 
is quadratic on the amount of time slices. We managed to optimize this part by reducing each iteration to 
about 10 double precision operations. 

In essence, this algorithm is not slower than standard binomial trees and it is very simple to 
implement. 

4. CALCULATION OF PROBABILITY 

4.1. Direct Calculation of Probability 

We can calculate the probability density function by first recursively computing the probabilities of 
reaching each node of the tree. This can be performed efficiently thanks to the Markov property. To 
compute the density function we need to rescale these probabilities by the distance to the neighboring 
nodes: the more spread the nodes are, the lower the density. We can estimate the probability density as 
follows: First we compute the probability of reaching each final node of the tree. We do this 
incrementally by first computing the probabilities of reaching nodes in time slice 1, then time slice 2 and 





I . 



(16) 



C{Si, t-At) = e-'-'^'lpCiSi^u t) + qC{Si.i, t)] . 



(17) 
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so forth. At time slice 0, we know that the middle node has probability 1 of being reached and all the 
others have probabihty 0. We compute the probability of reaching a node as a sum of two contributions 
from the previous time slice. We reach the node with transition pu from the node below at the previous 
slice, and with transition pd from the node above. Each contribution is the product of the probability at 
the previous node times the transition to the current node. This formula is just a discretized version of the 
Chapman-Kolmogorov equation 

p(Xj, ti) = p(Xj_i,ti_i)puj_-i + pixj+i,ti_i)pdj+i . (18) 

Now that we have computed the absolute probabilities at the final nodes, we can give a proper 
prepoint-discretized estimation of the density by scaling the probabilities by the spread of the S values. 
For the upper half of the tree we divide the probability of each final node by the size of the lower adjacent 
interval in the mesh: density^ = PiK^i ~ ^1-2)- (Note: We use index 5',_2 because the binomial tree is 
constructed over a trinomial tree. In this way we can keep in memory all the nodes but only half of the 
nodes though are true final nodes.) If there is a final middle node we divide its probability by the average 
of sizes of the two adjacent intervals, that is: density, = /',/((5',+2 - 5,_2)/2). For the lower half of the 
mesh we divide the probability by the upper adjacent gap in the mesh: density,- = p,/(5',+2 - 5,). 

4.2. Numerical Derivatives of Expectation of Probability 

The probability P can be calculated as a numerical derivative with respect to strike Z of a European 
Call option, taking the risk-free rate r to be zero, given an underlying evaluated at time t = 0, with 
strike X, and other variables such as volatiUty a, cost of carry b, and time to expiration T suppressed here 
for clarity, C(So, 0; X), 

P[SiT)\Sito)] = P[X\Sito)] = l^ (19) 

S{T)=X dX^ 

This calculation of the probability distribution is dependent on the same conditions necessary for 
any tree algorithm, i.e., that enough nodes are processed to ensure that the resultant evaluations are a good 
representation of the corresponding Fokker-Planck equation, addressed above, and that the number of 
iterations in PATHTREE are sufficient for convergence. 

4.2.1. Alternative First Derivative Calculation of Probability 

An alternative method of calculating the probability P a a first-order numerical derivative, instead 
of as second-order derivative, with respect to X is to define a function Ch using the Heaviside step- 
function H{S, X) = I if S >X and otherwise, instead of the Max function at the time to expiration. 
This yields 

P[S{Tmto)] = P[X\S{to)] = - ^ (20) 

S{T)=X uJL 

Sometimes this is numerically useful for sharply peaked distributions at the time of expiration, but we 
have found the second derivative algorithm above to work fine with a sufficient number of epochs. 

Our tests verify that the three methods above give the same density. We consider the numerical- 
derivative calculations a very necessary baseline to determine the number of epochs required to get 
reasonable accuracy. 

4.2.2. Oscillatory Corrections 

Fig. 1 illustrates the importance of including oscillatory corrections in any binomial tree algorithm. 
When these are included, it is easy to see the good agreement of the BS PATHTREE and OU PATHTREE 
models. 

4.3. Comparison to Exact Solutions 

Fig. 2 gives the calculated probabihty distribution for the BS and OU models, compared to their 
exact analytic solutions. 



Probability tree . 



-8- 



Ingber, Chen, Mondescu, Muzzall, Renedo 



5. CONCLUSION 

We have developed a path-integral based binomial PATHTREE algorithm that can be used in a 
variety of stochastic models. This algorithm is simple, fast and can be applied to diffusion processes with 
quite arbitrarily nonlinear drifts and diffusions. 

As expected, this PATHTREE algorithm is not as strong as PATHINT [4], as PATHINT can include 
details of an extremely high dimensional tree with complex boundary conditions. 

For PATHINT, the time and space variables are determined independently. I.e., the ranges of the 
space variables are best determined by first determining the reasonable spread of the distribution at the 
final time epoch. For PATHTREE just one parameter, the number of epochs A^, determines the mesh for 
both time and the space variables. This typically leads to a growth of the tree, proportional to V^, much 
faster than the spread of the distribution, so that much of the calculation is not relevant. 

However, this PATHTREE algorithm is surprisingly robust and accurate. Similar to PATHINT, we 
expect its accuracy to be best for moderate-noise systems. 
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FIGURE CAPTIONS 

FIG. 1. The oscillatory correction, an average of and A'^ + 1 iteration solutions, provides a simple 
and effective fix to the well-known oscillations inherent to binomial trees. The uncorrected Black-Scholes 
binomial tree (a) can be compared to the Black-Scholes tree with oscillatory correction (b). In (c), the 
Omstein-Uhlenbeck binomial tree also be robustly corrected as shown in (d). The BS PATHTREE model 
shown in (e) can be compared to the Black-Scholes case shown in (b). The OU PATHTREE model (f) is 
equivalent to the Omstein-Uhlenbeck model in (d). Parameters used in these calculations are: S = 50.0, X 
= 55.0, r = 1.0, r = 0.0675, b^O,a^ 0.20, and = 300. 

FIG. 2. Probability distributions for the PATHTREE binomial model as described in the text. In 
(a), bar graphs indicate OU PATHTREE agrees well with the exact Omstein-Uhlenbeck distribution 
shown in the black line. In (b), the bar graphs indicate BS PATHTREE agrees well with the exact Black- 
Scholes distribution shown in the black line. Parameters are the same as in Fig. I . 
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